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Results from direct numerical simulations for three dimensional Rayleigh-Benard con- 
vection in a cylindrical cell of aspect ratio 1/2 and Pr = 0.7 are presented. They span 
five decades of Ra from 2 x 10^ to 2 x 10^^. Good numerical resolution with grid spacing 
Kolmogorov scale turns out to be crucial to accurately calculate the Nusselt number, 
which is in good agreement with the experimental data by Niemela et ai, Nature, 404, 
837 (2000). In underresolved simulations the hot (cold) plumes travel further from the 
bottom (top) plate than in the fully resolved case, because the thermal dissipation close 
to the sidewall (where the grid cells are largest) is insufficient. We compared the fully 
resolved thermal boundary layer profile with the Prandtl-Blasius profile. We find that 
the boundary layer profile is closer to the Prandtl Blasius profile at the cylinder axis 
than close to the sidewall, due to rising plumes in that region. 



1. Introduction 

Turb ulent Rayleigh-Bena r d convection ( RBC) , continues to be a topic of intense re- 
search ( Ahlers et al. ( 2009h : Lohse fc Xial ( 2010l )). The system is relevant to numerous 
astro- and geo-physical phenomena, including convection in the arctic ocean, in Earth's 
outer core, in the interior of gaseous giant planets, and in the outer layer of the Sun. Thus 
the problem is of interest in a wide range of sciences, including geology, oceanography, 
climatology, and astrophysics. 

For given aspect ratio T = D/L [D is the cell diameter and L its height) and given 
geometry, the nature of RBC is determined by the Rayleigh number Ra ~ jSg/S.L^ / {kv) 
and by the Prandtl number Pr — v j k. Here, /? is the thermal expansion coefficient, g 
the gravitational acceleration, A = — Tj the difference between the imposed temper- 
atures Tb and Tt at the bottom and the top of the sample, respectively, and v and k the 
kinematic viscosity and the thermal diffusivity, respectively. When Pr is fixed, Ra is the 
only control parameter. There is still no universally ac cepted theory wha t the asymptoti c 
high Rayleigh numbe r Nu(Ra) relationship shoul d be (Kraichnaiil ( 1962 ): Spiegell (1971); 
Castaing all (|l989[ ): IShraima n fc Sig gial (Il990h: I Ahlers tt al\ (|2009[V ). and t he experi- 
menta l results are controversial [H e slot ei al. (Il987 ):IChavanne et aZj(|l997): Roche et al 



■20021) :lNiemela aZ.I (tioOOtbOOlh ^ lNiemela fc SreenivasanI (|2003f) : lNikolaenko et al\^20m : 
FunfschiUing et al\ (|2005l . l2009l )). 

For more mo d erate Ra up to 2 x 10^'^ previous direct numerical simulations (DNS) 



bv lAmati et al. ( 2005 ) in a three dimensional cylindrical cell of aspect ratio 1/2 with 
Pr = 0.7 showed a higher Nusselt Nu number than mea sured in experiments, see figure 
[2 To explain this discrepancy it was then suggested by IVerzicco fc Sreenivasan ( 2008f ) 
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Figure 1 . (a) Compensate d Nusselt number vs the Rayleigh number for Pr = 0.7. Purple stars 
are ex perimental data from|Nicmcla et aV. (2000) and the green squ ares are from | Chavanne et al\ 
l|200ll '). The DNS results from Verzicco & Camussi (2003) and A mati et al\ (|2005l ) are indicated 
in red and the present DNS results with the highest resolution are indicated by the black dots 
(the error is smaller than the dot size). The results of the underresolved simulations of this study 
are indicated by the blue dots, (b) Sketch of the grid geometry. The cells close to the sidewall 
are largest and therefore this region is least resolved. 



that the experimental conditions are closer to fixed flux condition s than fixed temper- 
ature conditions. However, recent two dimensional simulations bv I Johnston fc Doering 
( 20091 ) showed that Nu obtained in simulations with constant temperature and constant 
heat flux are identical when Ra > 5 x 10^. In this paper we show that the Nusselt num- 
ber obtained in the three dimensional simulations with constant temperature conditions 
is in good agreement with the experimental data, see flgure [TJ when the resolution is 
sufliciently high. 



2. Numerical method and results on the Nusselt number 

We numerically solved the three-dimensional Navier-Stokes equations within the Boussi- 
nesq approximation, 



Dm f Pr 

= -VP+ — 

Dt \Ra 



1/2 



DO _ 1 
In ^ {PrRay/^ 



(2.1) 
(2.2) 



with V • u = 0. Here z is the unit vector pointing in the opposite direction to gravity, 
D/ Dt = dt + M-y the material derivative, u the velocity vector (with no-slip boundary 
conditions at all walls), and the non-dimensional temperature, < 6* < 1. The equations 
have been made non-dimensional by using, next to L a nd A, the free-fall velocity U = 
\/ (5q^L. The numer i cal scheme is described in detail in Verzicco fc Orlandi ( 19961 ) and 
Verzicco fc Camussil (Il999i |2Q03[) . 

An important criterion in DNS simulations is that all the relevant flow scales, i.e. the 
Kolmogorov length r\ and the Batchelor length ry^, are properly resolved. Si nce Pr < 1, 
the smallest of these is r\ and we d etermined rjhy rj/L k. Tr{Pr'^ /RaNu)^^^ (Grotzbach 
(|l983l ): IVerzicco fc Camussil (l2003l) 1. We used different grids to test the influence of the 
grid scales. In table [1] the largest grid scale imax is compared to the Kolmogorov scale 
r] for each simulation. Note that their ratio is based on the global criterion, assuming a 



Radial boundary layer structure and Nusselt number in Rayleigh-Benard convection 3 



uniform distribution of the dissipation rates, in contrast to the observed peaking of the 
dissipation rates close to the (side) wall (see figure [2]). This means that the resolution in 
the bulk is better than the indicated value, however worse close to the sidewalk The grid 
density near the plates has been enhanced to keep a sufficient number of nodes in the 
thermal boundary layer (BL) where the dissipation rates are high, see the column " Ae" in 
table [H We calculate Nu as volume average and also by using the temperature gradients 
at the bottom and top plate. The volume average is calc ulated from the defi n ition of 
the Nusselt number Nu = [{u^e)A - ^83(9) a)/ kAL'^ (jVerzicco fc Camussil (|l999l) ). 
In addition, we average over the entire volume and time. The averages of the three 
methods, i.e. the volume average and the averages based on the temperature gradients 
at the bottom and top plate, are determined over at least 400 dimensionless time units. 
The value Nu in table [T] gives the average value of these three. We also determined Nu 
over the last half of our simulations, see the column " Nuh" in table [T] These values are 
within 1% of the value determined over the whole simulation, showing that our results 
arc well converged. The maximum difference in Nu obtained from the three methods, i.e. 
volume average and using the temperature gradients at the plates, is given in the column 
"conv" in table [TJ We simulated 200 dimensionless time units before we started to collect 
data to be sure to have reached the statistically stationary state. This is confirmed by 
comparing the statistics of the last part of the simulation used for the initialization and 
the part of the simulation used for the actual data collecting. In simulations where the 
field obtained at a lower Ra (or a new random field) is used as initial condition, we 
observe a small overshoot in Nu, before it settles to it statistically stationary value, 
The long initialization runs are thus used to eliminate this effect. This is double checked 
by the convergence of the three different methods to calculate Nu. We never noticed a 
(significant) drift in Nu when the results of the three methods arc within 1% of each 
other. For the four most demanding simulations, i.e. the bottom four cases in table [TJ 
the criteria for time averaging had to be relaxed due to the limited CPU time available. 
Therefore we averaged these cases for 100 dimensionless time units (300 time units for 
the simulation at Ra = 2 x 10^" on the 385 x 257 x 1025 grid). Wc note that the error bars 
given in figure [T] show the error in Nusselt obtained on that particular grid. The effect of 
the grid resolution on the Nusselt results for the highest Ra numbers is not estimated. 

Since most simulations are started from a interpolated field obtained at a lower Ra, we 
recomputed Nu for Ra = 2 x 10*^ on the 193 x 65 x 257 grid with a new flow field to rule 
out the effect of hysteresis on the obtained Nusselt results. The result is shown in italics 
in table [1] and is in excellent agreement with the original result. This confirms that our 
initialization and data collecting runs are long enough to eliminate the hysteresis effect 
on the Nusselt results. The simulations at Ra = 2 x 10^° have different initial conditions, 
i.e. different flow fields obtained at lower Ra are used as initial condition, and here we also 
observe a good agreement considering the different grids that are used. Figure [1] shows 
that the DNS data converge to the experimental data when the resolution is increased. 
To obtain accurate results for the Nusselt number in DNS simulations the grid spacing 
has to < 77 in the whole domain. We note that in simulations for Pr > 1 it is important 
to properly resolve the Batchelor scale rjT over the whole domain, since tjt < r] when 
Pr > 1 . 



3. Dissipation rates, temperature distribution functions, and 
boundary layers 

Another way to calculate Nu is to look at the two exact global relations for the volume 
averaged kinetic and thermal energy dissipation rates (e„) = v^{Nu — l)RaPr^^/L'^, 
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Table 1 . The columns Irom left to right indicate Ra; the number of grid points in the azimuthal, 
radial, and axial direction {Nq x Nr x N^), the Nusselt number (A'^m) obtained after averaging the 
results of the three methods (see text) using the whole simulation length, the Nusselt number 
(Nuh) after averaging the results of the three methods using the last half of the simulation, 
the maximum difference between the three methods (Conv), the number (Nbl) of points in 
the thermal BL, the maximum grid scale compared to the Kolmogorov scale estimated by the 
global criterion {Imax/rj). The last two columns give the Nusselt number derived from the volume 
averaged kinetic (eu) and thermal (eg) dissipation rates compared to Nu indicated in column 
three. The italic line indicates a simulation started with a new flow fleld. 



Ra 


Ng X Nr X Nz 


Nu 


Nuh 


Conv 


Nbl 


^max 1 ^ 










Nu 


JVu 


2 X 10*^ 


97 X 49 X 129 


10.85 


10.92 


0.32 


% 


18 


0.42 






2 X 10^ 


129 X 49 X 193 


20.52 


20.56 


0.36 


% 


17 


0.66 






2 X 10* 


97 X 49 X 193 


40.57 


40.71 


0.02 


% 


10 


1.84 


1.0062 


0.8528 


2 X 10* 


193 X 65 X 257 


39.42 


39.52 


0.02 


% 


13 


0.92 


0.9901 


0.9001 


2 X 10* 


257 X 97 X 385 


39.41 


39.10 


0.79 


% 


19 


0.70 


0.9924 


0.9411 


2 X lO'' 


129 X 65 X 257 


89.07 


88.25 


0.02 


% 


6 


3.01 


1.0000 


0.7317 


2 X l(f 


193 X 65 X 257 


84.49 


84.46 


0.45 


% 


7 


1.99 


1.0009 


0.7638 


2 X 10^ 


193 X 65 X 257 


84.10 


83.66 


0.51 


% 


7 


1.98 


0.9993 


0.7625 


2 X 10'' 


385 X 97 X 385 


79.75 


78.70 


0.70 


% 


10 


1.15 


0.9974 


0.8693 


2 X 10^° 


129 X 97 X 385 


201.08 


201.21 


1.01 


% 


12 


6.56 


1.0058 


0.8661 


2 X 10"^ 


513 X 129 X 513 


171.79 


169.58 


2.09 


% 


19 


1.59 


0.9981 


0.9234 


2 X 10"^ 


385 X 257 X 1025 


173.13 


173.30 


0.98 


% 


29 


2.12 






2 X 10" 


769 X 193 X 769 


387.07 


387.53 


2.18 


% 


16 


2.31 






2 X 10" 


769 X 257 X 1025 


373.64 


368.88 


2.03 


% 


18 


2.28 


0.9828 


0.8910 



and (eg) = kA'^Nu/L^, respectively ( Shraiman fc Siggia ( 1990[ )). We have calculated the 
azimuthally and time averaged energy dissipation rate eui'st) = z^jVup and the thermal 
dissipation rate e6i(~af ) = k| V0p. Figure[2]compares the difference between the dissipation 
rates obtained in the fully resolved and the underresolved simulations and reveals a 
higher thermal dissipation rate for the fully resolved simulations as it is calculated from 
the (temperature) gradients. In the underresolved simulations the gradients are smeared 
out and therefore and eg are underestimated. To check the resolution, we calculated 
Cu and eg from the respective gradients and compared it with the values obtained from 
above global exact relations. Tabic [T] shows that for e„ the relation is basically satisfied, 
whereas for eg the difference is considerable when the simulation is underresolved, and 
even for the best resolved cases there is 6 — 8% too little dissipation. Testing above exact 
relations seems to be the best way to verify the grid resolution. 

The vertical heat flux concentrates i n the plume-domina ted sidewall region where the 
vertical velocity reaches its maximum (jShang et al. I (I2OO8I) ). Therefore it is very impor- 
tant to properly resolve the region close to the sidewall. However, figure [2] reveals that in 
the underresolved simulations the region close to the sidewall is least resolved (red areas 
in the plot where the thermal dissipation rates are compared (right plot)), as there the 
grid boxes are largest, due to the cylindrical geometry of the grid (see figure [Ud). When 
the resolution is insufficient close to the sidewall, the plumes in this region, important 
for the heat transfer, are not properly resolved and not sufficiently dissipated. Therefore 
too much heat reaches the other side and correspondingly Nu is overestimated in these 
underresolved simulations. Supplementary movies reveal the dynamics of the system for 
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Figure 2. Dimensionless kinetic (upper) and thermal (lower) dissipation rates at Ra — 2 x 10^. 
The upper row gives = euL^/U'^ and the lower row ee — eg [//(A^i). The left column 
indicates the dimensionless kinetic and thermal eg dissipation rates for the high resolution case 
(385 X 97 X 385). The middle column gives — (upper plot), and ef^ — (lower plot), where 
the superscripts H and L, respectively, mean the data obtained from the high (385 x 97 x 385) 
and low resolution simulations (129 x 65 x 257). The rightmost column gives (e^ —^u)/^u (upper 
plot) and (ef — 'eg)/(:f (lower plot). The difference for the thermal dissipation rates between 
the fully resolved and the underresolved simulations is largest (in absolute values) close to the 
sidewalk 



the different grid resolutions. Movie 1 shows the temperature field close to bottom plate 
and movie 2 the temperature field at mid height. Note that the smoothness of the under- 
resolved simulations is insufficient to capture all characteristics of the flow represented 
in the high resolution simulation. 

To further investig ate the influ ence of the grid resol ution, we calc ulated azimut hally av- 

eraged PDF s (see also|Emran fc S chumacher (2008.): .Kunnen et all ((2008, ): .Shishkina fc Wagner 
( 2007ll2008h : iKaczorowski fc Wagneii (|2009h ) of the temperature averaged over 3000 di- 
mensionless time units for Ra = 2 x 10®, comparing the underresolved case (97 x 49 x 193) 
with the fully resolved one (193 x 65 x 257). Figure [3] shows that the temperature PDFs 
at mid height and at a distance Xf (thermal BL based on the slope) from the plates 
have longer tails in the underresolved simulation than in the fully resolved one. Again 
the reason lies in the rising (falling) plumes from the bottom (top) plate which are not 
properly dissipated in the underresolved simulations and therefore travel further from the 
plates. The comparison with the PDF obtained using half of the time series reveals that 
the differences in the PDFs are not due to a lack of averaging, but due to insufficient grid 
resolutions. We note that we observe similar differences at other radial positions, only 
the averaging around the cylinder axis (r = 0) leads to not fully converged results due to 
the geometry. In figure|3]we show the effect of the grid resolution on the flatness obtained 
at mid height for fully resolved and underresolved simulations. Comparison between the 
solid and dashed lines show that the data are converged close to the sidewall where the 
statistics is best due to geometric reason. Comparison between black (well resolved) and 
red (underresolved) reveals that the insufficiently dissipated plumes mainly close to the 
sidewall leads to too large flatnesses in the underresolved simulations. 

Although the bulk is turbulent, scalingwise the BLs still behave in a laminar way 
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(c) (d) 

Figure 3. a) A sketch showing the locations (crosses) of the azimuthally averaged temperature 
PDFs, shown in figure |3]b, c, and d, for Ra = 2 x 10* obtained on different grids. The radial 
position is 0.2342L for the underresolved (97 x 49 x 193) and 0.2314L for the fully resolved 
(193 X 65 X 257) simulations. The temperature PDF for the fully resolved simulations averaged 
over 3000 dimensionless time units is indicated in black. The green line indicates the result using 
half of the time series. The temperature PDF averaged over 3000 dimensionless time units for 
the underresolved simulations is indicated in blue, and the red indicates the result using half of 
the time series, b) Temperature PDF at mid height, c) Temperature PDF at the distance Xf 
from the bottom plate, d) Temperature PDF at the distance Xf from the top plate. 



due to the small BL Reynolds number ( Ahlers et al. ( 20091 )). Therefore we compare the 



thermal B L profile obtained from the simulations with the Prandtl-Blasius (PB) profile, 
as done bv lSugivama et 'd\ (l2009l) for 2D RB simulations. The temperature gradient of 
the PB profile is matched to the temperature gradient obtained in the high resolution 
simulation. The temperature profile obtained in the simulations best matches the PB 
profile around the cylinder axis (r = 0). Close to the sidewall the agreement is worse due 
to the rising (falling) plumes in this region. We now compare the difference between the 
PB profile and the result obtained from the simulation for different Ra. We determine, 
at the cylinder axis, {9 sim-S p b) / {^S p b) for the bottom BL and {0 p b-9 sim) / {0 p b) for 
the top BL. Here Osim is the mean temperature at a distance \f from the plate and 
OpB the temperature according to PB at this height, after having matched the gradient 
at the plate to the simulation data. If the simulation exactly matched PB (e.g. for very 
small Ra), this expression would be zero. In contrast, here for Ra = 2 x 10^ (2 x 10^, 
2 X 10^°), it is 0.103, (0.130, 0.149). As expected, the expression is smaller for the lower 



Ra numbers. We perform the same procedure for our previous results of IZhong et 
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Figure 4. Flatness of the temperature PDF at mid height for Ra = 2 x 10* for the underresolved 
(red, 97 X 49 X 193) and the fully resolved simulations (black, 193 x 65 x 257). The sohd lines 
indicate the result after averaging over 3000 dimensionless time units and the dashed lines the 
result after averaging over 1600 dimensionless time units. Both simulations are started from the 
same initial field obtained at a lower Ra and the data collecting is started when each simulation 
has reached the statistically stationary state. 




Figure 5. Azimuthally averaged temperature profiles obtained from the simulations at different 
grids, a) i?a = 2 X 10^ for the grids 97 x 49 x 193 (red), 193 x 65 x 257 (blue), 257 x 97 x 385 
(black) and b) iia = 2 x lO'-* for the grids 129 x 65 x 257 (red), 193 x 65 x 257 (blue), 385 x 97 x 385 
(black). The diamonds indicate the axial position of the grid points. The solid lines show the 
temperature profile at the cylinder axis (r — 0) and the dashed lines at the radial position 
0.2251/. The green line indicates the PB profile matched to the temperature gradient at the 
cylinder axis (r = 0) of the high resolution simulation. The insets show the temperature profile 
from the highest resolution data over a larger axial range. Here the solid line indicates the profile 
at the axis and the dashed line the temperature profile at the radial position 0.2251/. 

( 20091 ) at i?a = 1 x 10** with T = 1 and now different Pr. For Pr = 0.7, Pr = 6.4, 
and Pr = 20 we now obtain 0.099, 0.040, and 0.033, respectively. Now the expression is 
closest to zero at the higher Pr, as then the Re number is lower and PB holds better. 

Figure [6] shows the radial dependence of Xf , A™" (thermal BL thickness based on 
maximum rms value). A™* (kinetic BL thickness based on maximum azimuthal rms 
velocity), and A^" (kinetic BL thickness defined as the axial position of the maximum 
kinetic energy dissipation rate, multiplied by 2). A™* is widely used in literature to define 
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the kinetic BL thickness; however, this definition overestimates the kinetic BL thickness. 
A^" defines the BL as the region where the kinetic dissipation is highest and it is this 
region where a particular good resolution is required. Such defined kinetic BL thickness 
now well agrees with that of the thermal BL, A^" w Xf, as from PB theory expected for 
the kinetic BL, once Pr ~ 1. Figure [6] shows that both BLs become thicker closer to the 
sidewalk This is due to the plumes traveling along the sidewall and lower velocities very 
close to the sidewall. Thus the enhanced grid resolution in the vertical direction near 
the plates is most important around the cylinder axis (r = 0). In contrast, the azimuthal 
(and radial) resolution arc most important when properly resolving the flow close to 
the sidewall. Note that the difference in the BL thicknesses between the fully resolved 
and underresolved simulations is largest close to the sidewall, demonstrating that this is 
indeed a delicate region from a resolution point of view. 



4. Conclusions 

In summary, the high-resolution results using constant temperature conditions arc in 
good agreement with the experimental data, sec figure [TJ It thus turns out that a good 
resolution of imax ^ Vt is crucial to obtain reliable results for Nusselt number. Close 
attention has to be given to the resolution used in all directions (azimuthal, radial, and 
axial). In particular, the azimuthal resolution is crucial to guarantee sufficient plume 
dissipation close to the sidewall. In underresolved simulations the exact relation eg = 
kA'^ Nu/ L"^ for the thermal dissipation rate does not hold. Hot (cold) plumes travel 
further from the bottom (top) plate than in the fully resolved simulation, because they 
are not sufficiently dissipated. We also showed that there is a strong radial dependence 
of the BL structures. At the cylinder axis (r=0) the temperature profile obtained in the 
simulations agrees well with the PB case, whereas close to the sidewall the agreement is 
worse due to rising (falling) plumes in this region. 

The effect of changing the constant temperature condition at the bottom plate to a 
constant heat flux condition will be discussed in detail in a forthcoming publication. 
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